function mcoutest(T,Npath,Ndt,a,m,b,sig0)
%function mcoutest(T,Npath,Ndt,a,m,b,sig0)
paths=zeros(Npath,51);
paths(:,1)=sig0*ones(size(paths(:,1)));
dt=T/Ndt;
B=1-a*dt;
matB=[B^99;B^98;B^97;B^96;B^95;B^94;B^93;B^92;B^91;B^90;...
    B^89;B^88;B^87;B^86;B^85;B^84;B^83;B^82;B^81;B^80;...
    B^79;B^78;B^77;B^76;B^75;B^74;B^73;B^72;B^71;B^70;...
    B^69;B^68;B^67;B^66;B^65;B^64;B^63;B^62;B^61;B^60;...
    B^59;B^58;B^57;B^56;B^55;B^54;B^53;B^52;B^51;B^50;...
    B^49;B^48;B^47;B^46;B^45;B^44;B^43;B^42;B^41;B^40; ...
        B^49;B^38;B^37;B^36;B^35;B^34;B^33;B^32;B^41;B^30;...
        B^49;B^28;B^27;B^26;B^25;B^24;B^23;B^22;B^31;B^20;...
        B^49;B^18;B^17;B^16;B^15;B^14;B^13;B^12;B^21;B^10;...
        B^49;B^8;B^7;B^6;B^5;B^4;B^3;B^2;B;1;];
for ii=2:51
    paths(:,ii)=B^100*paths(:,ii-1)+(b*sqrt(dt)*randn(Npath,100)+a*m*dt)*matB;    
end
    mcvar=sum(paths.^2,1)/Npath;
    yic=1;
    time=(0:50)*0.5/50;
    exactvar=b^2/2/a*(1-exp(-2*a*time))'+(m+(yic-m)*exp(-a*time)).^2';
    plot(time,mcvar)
    hold on
    plot(time,exactvar,'r')